# The politics of integrating health systems and public programs
# N.M. Kavanagh, A. Menon, A. McIntyre
# Code primarily authored by A. Menon

# Last updated: January 13, 2026

# Load packages
library(here)         # Working directory
library(tidyverse)    # Analysis tools
library(dplyr)        # Analysis tools
library(survey)       # Modeling tools
library(ggplot2)      # Graphing tools
library(ggpubr)       # Graphing tools
library(ggtext)       # Graphing tools
library(patchwork)    # Graphing tools

##############################################################################
# Clean ESS dataset
##############################################################################

# Read in ESS dataset
ess <- read.csv("ESS download/ESS waves 1-11 dataset.csv")

# Trust in politicians
ess <- ess %>% mutate(
  trstplt_d = case_when(
    trstplt %in% c(0:10) ~ trstplt))

# Trust in parliament
ess <- ess %>% mutate(
  trstprl_d = case_when(
    trstprl %in% c(0:10) ~ trstprl))

# Trust in political parties
ess <- ess %>% mutate(
  trstprt_d = case_when(
    trstprt %in% c(0:10) ~ trstprt))

# Trust in European Parliament
ess <- ess %>% mutate(
  trstep_d = case_when(
    trstep %in% c(0:10) ~ trstep))

# Trust in legal system
ess <- ess %>% mutate(
  trstlgl_d = case_when(
    trstlgl %in% c(0:10) ~ trstlgl))

# Trust in police
ess <- ess %>% mutate(
  trstplc_d = case_when(
    trstplc %in% c(0:10) ~ trstplc))

# Trust in United Nations
ess <- ess %>% mutate(
  trstun_d = case_when(
    trstun %in% c(0:10) ~ trstun))

# Satisfaction with how democracy works
ess <- ess %>% mutate(
  stfdem_d = case_when(
    stfdem %in% c(0:10) ~ stfdem))

# Satisfaction with present economy
ess <- ess %>% mutate(
  stfeco_d = case_when(
    stfeco %in% c(0:10) ~ stfeco))

# Satisfaction with national government
ess <- ess %>% mutate(
  stfgov_d = case_when(
    stfgov %in% c(0:10) ~ stfgov))

# State of health services
ess <- ess %>% mutate(
  stfhlth_d = case_when(
    stfhlth %in% c(0:10) ~ stfhlth))

# State of education system
ess <- ess %>% mutate(
  stfedu_d = case_when(
    stfedu %in% c(0:10) ~ stfedu))

# Generate analysis weights
# Combine post-stratification and population weights
ess$analysis_wt <- ess$pspwght * ess$pweight

# Dummy variable necessary for survey table
ess$one <- 1

##############################################################################
# Main graph of descriptive statistics (based on listwise deletion)
##############################################################################

# Get Ns for respondents with complete data
ess %>%
  filter_at(vars(trstplt_d, trstprl_d, trstprt_d, trstep_d,
                 trstlgl_d, trstplc_d, trstun_d, stfdem_d,
                 stfhlth_d, stfeco_d, stfedu_d, stfgov_d),
            all_vars(!is.na(.))) %>%
  summarise(n_resp  = n(),
            n_cntry = length(unique(cntry)),
            waves   = unique(essround))

# Exclude wave 1 since not all outcomes assessed
ess <- ess %>% filter(essround != 1)

# Generate survey design
design <- svydesign(id=~1, weights=~analysis_wt, data=ess)

# Generate weighted means for each outcome
# Note: This procedure enforces listwise deletion of missing data
sum_stats <- svyby(~ trstplt_d + trstprl_d + trstprt_d + trstep_d +
                     trstlgl_d + trstplc_d + trstun_d + stfdem_d +
                     stfhlth_d + stfeco_d + stfedu_d + stfgov_d, ~one,
                   design, svymean, na.rm=TRUE, vartype="ci")

# Restructure table
ss_long <- sum_stats %>%
  
  # Pivot to long form
  pivot_longer(cols = everything(), names_to = "variable", values_to = "value") %>%
  
  # Rename variables
  mutate(
    type = case_when(
      grepl("^ci_l\\.", variable) ~ "ci_low",
      grepl("^ci_u\\.", variable) ~ "ci_high",
      TRUE ~ "mean"),
    var_clean = gsub("^ci_[lu]\\.", "", variable)) %>%
  
  # Select desired variables
  select(var_clean, type, value) %>%
  
  # Revert to wide format
  pivot_wider(names_from = type, values_from = value) %>%
  
  # Remove dummy variable
  filter(var_clean != "one") %>%
  
  # Reorder by mean
  arrange(desc(mean))

# Labels of variables
labels <- c(
  trstplt_d = "Trust in politicians",
  trstprl_d = "Trust in parliament",
  trstprt_d = "Trust in political parties",
  trstep_d  = "Trust in European Parliament",
  trstlgl_d = "Trust in legal system",
  trstplc_d = "Trust in police",
  trstun_d  = "Trust in United Nations",
  stfdem_d  = "Satisfaction with democracy",
  stfeco_d  = "Satisfaction with economy",
  stfgov_d  = "Satisfaction with government",
  stfhlth_d = "State of health services",
  stfedu_d  = "State of education system"
  )

# Add variable labels
# Include highlight indicator for health services
ss_long <- ss_long %>%
  mutate(label     = labels[var_clean],
         highlight = ifelse(var_clean == "stfhlth_d", TRUE, FALSE))

# Subset to variables with each scale
ss_state        <- ss_long[grepl("State", ss_long$label), ]
ss_trust        <- ss_long[grepl("Trust", ss_long$label), ]
ss_satisfaction <- ss_long[grepl("Satisfaction", ss_long$label), ]  

# Panel 1: State of various institutions
p_state <- ggplot(ss_state, aes(x = reorder(label, mean), y = mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2) +
  geom_text(aes(label = format(round(mean, 2), 2)),
            hjust = -0.8, size = 3) +
  coord_flip(clip = "off") +
  theme_minimal() +
  theme(text         = element_text(face = "bold"),
        axis.title.y = element_text(face="bold"),
        plot.margin  = margin(10, 70, 10, 20)) +
  labs(x = "", y = "\nMean reported state of various institutions") +
  scale_y_continuous(
    limits = c(0, 10),
    breaks = c(0,2,4,6,8,10),
    labels = c("0\nExtremely\nbad","2","4","6","8","10\nExtremely\ngood"))

# Add highlight for health services
p_state <- p_state +
  geom_point(
    data  = subset(ss_state, highlight),
    shape = 22, size = 9, fill = NA, color = "red", stroke = 1.5)

# Panel 2: Satisfaction with various institutions
p_satisfaction <- ggplot(ss_satisfaction, aes(x = reorder(label, mean), y = mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2) +
  geom_text(aes(label = format(round(mean, 2), 2)),
            hjust = -0.8, size = 3) +
  coord_flip(clip = "off") +
  theme_minimal() +
  theme(text         = element_text(face = "bold"),
        axis.title.y = element_text(face = "bold"),
        plot.margin  = margin(10, 70, 10, 20)) +
  labs(x = "", y = "\nMean reported satisfaction with various institutions") +
  scale_y_continuous(
    limits = c(0, 10),
    breaks = c(0,2,4,6,8,10),
    labels = c("0\nExtremely\ndissatisfied","2","4","6","8","10\nExtremely\nsatisfied"))

# Panel 3: Trust in various institutions
p_trust <- ggplot(ss_trust, aes(x = reorder(label, mean), y = mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2) +
  geom_text(aes(label = format(round(mean, 2), 2)),
            hjust = -0.8, size = 3) +
  coord_flip(clip = "off") +
  theme_minimal() +
  theme(text         = element_text(face = "bold"),
        axis.title.y = element_text(face = "bold"),
        plot.margin  = margin(10, 70, 10, 20)) +
  labs(x = "", y = "\nMean reported trust in various institutions") +
  scale_y_continuous(limits = c(0, 10),
                     breaks = c(0,2,4,6,8,10),
                     labels = c("0\nNo trust\nat all","2","4","6","8","10\nComplete\ntrust"))

# Combine into single panel
g <- p_state / p_satisfaction / p_trust + plot_layout(heights = c(0.6, 0.9, 2))

# Export figure
ggsave(plot=g, file="trust_statisfaction_all_main.tif",
       width=6, height=6, units='in', dpi=300)

##############################################################################
# Supplemental graph based on pairwise deletion
##############################################################################

# Get point estimates with pairwise deletion of missing data
sum_stats_pair <- ess %>%
  summarise(trstplt_d = weighted.mean(trstplt_d, w=analysis_wt, na.rm=T),
            trstprl_d = weighted.mean(trstprl_d, w=analysis_wt, na.rm=T),
            trstprt_d = weighted.mean(trstprt_d, w=analysis_wt, na.rm=T),
            trstep_d  = weighted.mean(trstep_d,  w=analysis_wt, na.rm=T),
            trstlgl_d = weighted.mean(trstlgl_d, w=analysis_wt, na.rm=T),
            trstplc_d = weighted.mean(trstplc_d, w=analysis_wt, na.rm=T),
            trstun_d  = weighted.mean(trstun_d,  w=analysis_wt, na.rm=T),
            stfdem_d  = weighted.mean(stfdem_d,  w=analysis_wt, na.rm=T),
            stfhlth_d = weighted.mean(stfhlth_d, w=analysis_wt, na.rm=T),
            stfeco_d  = weighted.mean(stfeco_d,  w=analysis_wt, na.rm=T),
            stfedu_d  = weighted.mean(stfedu_d,  w=analysis_wt, na.rm=T),
            stfgov_d  = weighted.mean(stfgov_d,  w=analysis_wt, na.rm=T))

# Get Ns by outcome
n_stats_pair <- ess %>%
  summarise(trstplt_d = sum(!is.na(trstplt_d)),
            trstprl_d = sum(!is.na(trstprl_d)),
            trstprt_d = sum(!is.na(trstprt_d)),
            trstep_d  = sum(!is.na(trstep_d)),
            trstlgl_d = sum(!is.na(trstlgl_d)),
            trstplc_d = sum(!is.na(trstplc_d)),
            trstun_d  = sum(!is.na(trstun_d)),
            stfdem_d  = sum(!is.na(stfdem_d)),
            stfhlth_d = sum(!is.na(stfhlth_d)),
            stfeco_d  = sum(!is.na(stfeco_d)),
            stfedu_d  = sum(!is.na(stfedu_d)),
            stfgov_d  = sum(!is.na(stfgov_d)))
min(n_stats_pair); max(n_stats_pair)

# Restructure table
ss_long_pair <- sum_stats_pair %>%
  
  # Pivot to long form
  pivot_longer(cols = everything(), names_to = "var_clean", values_to = "mean") %>%
  
  # Reorder by mean
  arrange(desc(mean))

# Add variable labels
# Include highlight indicator for health services
ss_long_pair <- ss_long_pair %>%
  mutate(label     = labels[var_clean],
         highlight = ifelse(var_clean == "stfhlth_d", TRUE, FALSE))

# Subset to variables with each scale
ss_state_pair        <- ss_long_pair[grepl("State", ss_long_pair$label), ]
ss_trust_pair        <- ss_long_pair[grepl("Trust", ss_long_pair$label), ]
ss_satisfaction_pair <- ss_long_pair[grepl("Satisfaction", ss_long_pair$label), ]  

# Panel 1: State of various institutions
p_state_pair <- ggplot(ss_state_pair, aes(x = reorder(label, mean), y = mean)) +
  geom_point() +
  geom_text(aes(label = format(round(mean, 2), 2)),
            hjust = -0.8, size = 3) +
  coord_flip(clip = "off") +
  theme_minimal() +
  theme(text         = element_text(face = "bold"),
        axis.title.y = element_text(face="bold"),
        plot.margin  = margin(10, 70, 10, 20)) +
  labs(x = "", y = "\nMean reported state of various institutions") +
  scale_y_continuous(
    limits = c(0, 10),
    breaks = c(0,2,4,6,8,10),
    labels = c("0\nExtremely\nbad","2","4","6","8","10\nExtremely\ngood"))

# Add highlight for health services
p_state_pair <- p_state_pair +
  geom_point(
    data  = subset(ss_state_pair, highlight),
    shape = 22, size = 9, fill = NA, color = "red", stroke = 1.5)

# Panel 2: Satisfaction with various institutions
p_satisfaction_pair <- ggplot(ss_satisfaction_pair, aes(x = reorder(label, mean), y = mean)) +
  geom_point() +
  geom_text(aes(label = format(round(mean, 2), 2)),
            hjust = -0.8, size = 3) +
  coord_flip(clip = "off") +
  theme_minimal() +
  theme(text         = element_text(face = "bold"),
        axis.title.y = element_text(face = "bold"),
        plot.margin  = margin(10, 70, 10, 20)) +
  labs(x = "", y = "\nMean reported satisfaction with various institutions") +
  scale_y_continuous(
    limits = c(0, 10),
    breaks = c(0,2,4,6,8,10),
    labels = c("0\nExtremely\ndissatisfied","2","4","6","8","10\nExtremely\nsatisfied"))

# Panel 3: Trust in various institutions
p_trust_pair <- ggplot(ss_trust_pair, aes(x = reorder(label, mean), y = mean)) +
  geom_point() +
  geom_text(aes(label = format(round(mean, 2), 2)),
            hjust = -0.8, size = 3) +
  coord_flip(clip = "off") +
  theme_minimal() +
  theme(text         = element_text(face = "bold"),
        axis.title.y = element_text(face = "bold"),
        plot.margin  = margin(10, 70, 10, 20)) +
  labs(x = "", y = "\nMean reported trust in various institutions") +
  scale_y_continuous(limits = c(0, 10),
                     breaks = c(0,2,4,6,8,10),
                     labels = c("0\nNo trust\nat all","2","4","6","8","10\nComplete\ntrust"))

# Combine into single panel
g_pair <- p_state_pair / p_satisfaction_pair / p_trust_pair + plot_layout(heights = c(0.6, 0.9, 2))

# Export figure
ggsave(plot=g_pair, file="trust_statisfaction_all_pairwise.tif",
       width=6, height=6, units='in', dpi=300)
